# Turn on the gmri font for plots
showtext::showtext_auto()

Gulf Of Maine Sea Surface Temperature Outlook

This report was created to track the sea surface temperature regimes for marine regions of interest to the Gulf of Maine Research Institute. The default region being a central snapshot of the Gulf of Maine.

Satellite sea surface temperature data used was obtained from the National Center for Environmental Information (NCEI). With all maps and figures displaying NOAA’s Optimum Interpolation Sea Surface Temperature Data.

Region Extent

The spatial extent for Gulf Of Maine is displayed below. This bounding box is the same bounding box coordinates used to clip the OISST data when constructing the time series data from the array.

# testing regions
# params <- list(); params$region <- "apershing_gulf_of_maine"

# File paths for various extents based on params$region
region_paths <- get_timeseries_paths(region_group = "gmri_sst_focal_areas")

# Load the bounding box for Andy's GOM to show they align
poly_path     <- region_paths[[params$region]][["shape_path"]]
region_extent <- st_read(poly_path, quiet = TRUE)


# Pull extents for the region for crop
crop_x <- st_bbox(region_extent)[c(1,3)]
crop_y <- st_bbox(region_extent)[c(2,4)]


# Zoom out for cpr extent
if(tolower(params$region) == "cpr gulf of maine"){
  crop_x <- c(-70.875, -65.375)
  crop_y <- c(40.375,   45.125)}


# one off shapes
# region_extent <- read_sf(paste0(res_path, ""))

# Full plot
ggplot() +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = gmri_cols("gmri blue"), 
          fill = gmri_cols("gmri blue"), alpha = 0.2, linetype = 2, size = 0.5) +
  map_theme +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = T)

Regional Timeseries

Area-specific time series are the most basic building block for relaying temporal trends. For any desired area (represented by a spatial polygon) we can generate a time series table of the mean sea surface temperature within that area for each day. Additionally, we can compare how observed temperatures correspond with the expected conditions based on a climatology using a specified reference period.

# Use {gmRi} instead to load timeseries to tye up loose ends
timeseries_path <- region_paths[[params$region]][["timeseries_path"]]
region_timeseries <- read_csv(timeseries_path, col_types = cols(), guess_max = 1e6)

# format dates
region_timeseries <- region_timeseries %>% 
  mutate(time = as.Date(time))

# Display Table of first 6 entries
tail(region_timeseries) %>% 
  mutate(across(where(is.numeric), round, 2)) %>% 
  rename(
    Date = time,
    `Sea Surface Temperature` = sst,
    `Day of Year` = modified_ordinal_day,
    `Climate Avg.` = sst_clim,
    `Climate SD` = clim_sd,
    `Temperature Anomaly` = sst_anom
  ) %>% gt() %>% 
    tab_header(
    title = md(paste0("**", tidy_name, " - Regional Sea Surface Temperature", "**")), 
    subtitle = paste("Temperature Unit: Celsius")) %>%
  tab_source_note(
    source_note = md("*Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.*") ) %>% 
  tab_source_note(md("*Climatology Reference Period: 1982-2011.*"))
Gulf Of Maine - Regional Sea Surface Temperature
Temperature Unit: Celsius
Date Sea Surface Temperature area_wtd_sst Day of Year Climate Avg. area_wtd_clim Climate SD Temperature Anomaly area_wtd_anom
2021-08-06 18.07 18.11 219 16.87 16.91 2.27 1.20 1.20
2021-08-07 18.54 18.57 220 16.90 16.93 2.26 1.64 1.64
2021-08-08 18.98 19.02 221 16.96 16.99 2.26 2.02 2.02
2021-08-09 19.04 19.08 222 17.06 17.10 2.24 1.98 1.98
2021-08-10 19.25 19.30 223 17.14 17.18 2.20 2.11 2.12
2021-08-11 19.48 19.53 224 17.20 17.23 2.21 2.28 2.29
Data Source: NOAA OISSTv2 Daily Sea Surface Temperature Data.
Climatology Reference Period: 1982-2011. 
# march 1st sst
mar1 <- region_timeseries %>% 
  filter(modified_ordinal_day == 61) %>% 
  distinct(sst_clim) %>% 
  pull(sst_clim)

Each of our Climatologies are currently set up to calculate daily averages on a modified year day, such that every March 1st and all days after fall on the same day, regardless of whether it is a leap year or not.

This preserves comparisons across calendar dates such-as: “The average temperature on march 1st is 4.2378197` for the reference period 1982 to 2011”

In these tables Sea Surface Temperature is the mean temperature observed for that date averaged across all cells within the area. Climate Avg. & Climate SD are the climate means and standard deviations for a 1982-2011 climatology. Temperature Anomaly is the daily observed sea surface temperature minus the climate mean.

Warming Rates

Regional warming trends below were calculated using all the available data for complete years beginning with 1982 through the end of 2020. The overlaid trend lines then track how warming has increased with time. A dotted line has been included to show how the global average temperature has changed during the same period.

Annual

# Summarize by year to return mean annual anomalies and variance
annual_summary <- region_timeseries %>% 
  mutate(year = year(time)) %>% 
  group_by(year) %>% 
  summarise(sst = mean(sst, na.rm = T),
            sst_anom = mean(sst_anom, na.rm = T), 
            .groups = "keep") %>% 
  ungroup() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))


# # Global Temperature Anomaly Rates
global_anoms <- read_csv(paste0(oisst_path, "global_timeseries/global_anoms_1982to2011.csv"))
global_anoms <- mutate(global_anoms, year = year(time))

global_summary <- global_anoms %>% 
  group_by(year) %>% 
  summarise(sst_anom = mean(sst, na.rm = T), .groups = "keep") %>% 
  ungroup() %>% 
  mutate(yr_as_dtime = as.Date(paste0(year, "-07-02")))
# Build Regression Equation Labels
lm_all <- lm(sst_anom ~ year, 
             data = filter(annual_summary, year %in% c(1982:2020))) %>% 
  coef() %>% 
  round(3)
lm_15  <- lm(sst_anom ~ year, 
             data = filter(annual_summary, year %in% c(2006:2020))) %>% 
  coef() %>% 
  round(3)
lm_global <- lm(sst_anom ~ year, 
                data = filter(global_summary, year %in% c(1982:2020))) %>% 
  coef() %>% 
  round(3)


# Convert yearly rate to decadal
decade_all    <- lm_all['year'] * 10
decade_15     <- lm_15['year'] * 10
decade_global <- lm_global["year"] * 10


# Equation to paste in
eq_all    <- paste0("y = ", lm_all['(Intercept)'], " + ", decade_all, " x")
eq_15     <- paste0("y = ", lm_15['(Intercept)'], " + ", decade_15, " x")
eq_global <- paste0("y = ", lm_15['(Intercept)'], " + ", decade_global, " x")
####  Annual Trend Plot  ####
ggplot(data = annual_summary, aes(yr_as_dtime, sst_anom)) +
  
  # Add daily data
  geom_line(data = region_timeseries,
            aes(time, sst_anom),
            alpha = 0.5, color = "gray") +
  
  # Overlay yearly means
  geom_line(alpha = 0.7) +
  geom_point(alpha = 0.7) +
  
  # Add regression lines 
  geom_smooth(data = filter(global_summary, year <= 2020),
            method = "lm", 
            aes(color = "1982-2020 Global Trend"),
            formula = y ~ x, se = F, 
            linetype = 3) +
  geom_smooth(data = filter(annual_summary, year <= 2020),
              method = "lm",
              aes(color = "1982-2020 Regional Trend"),
              formula = y ~ x, se = F,
              linetype = 2) +
  geom_smooth(data = filter(annual_summary, year %in% c(2006:2020)),
              method = "lm", 
              aes(color = "2006-2020 Regional Trend"),
              formula = y ~ x, se = F, 
              linetype = 2) +

  
  # Manually add equations so they show yearly not daily coeff
  geom_text(data = data.frame(), 
            aes(label = eq_all, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 2, color = gmri_cols("gmri blue")) +
  geom_text(data = data.frame(), 
            aes(label = eq_15, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 3.5, color = gmri_cols("orange")) +
  geom_text(data = data.frame(), 
            aes(label = eq_global, x = min(region_timeseries$time), y = Inf),
            hjust = 0, vjust = 5, color = gmri_cols("green")) +

  # Colors
  scale_color_manual(values = c(
    "1982-2020 Regional Trend" = as.character(gmri_cols("gmri blue")),
    "2006-2020 Regional Trend" = as.character(gmri_cols("orange")),
    "1982-2020 Global Trend"   = as.character(gmri_cols("green")))) +
  
  # labels + theme
  labs(x = "", 
       y = expression("Sea Surface Temperature Anomaly"~degree~C),
       caption = paste0("Regression coefficients reflect decadal change in sea surface temperature.
                         Anomalies calculated using 1982-2011 reference period.")) +
   theme(legend.title = element_blank(),
        legend.position = c(0.85, 0.1),
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent"),
        panel.grid = element_blank())

Quarterly

# Repeat for the seasons (equal quarters here*)
quarter_summary <- region_timeseries %>% 
  mutate(year = year(time),
         season = factor(quarter(time, fiscal_start = 1)),
         season = fct_recode(season, 
                             c("Jan 1 - March 31" = "1"), 
                             c("Apr 1 - Jun 30" = "2"), 
                             c("Jul 1 - Sep 30" = "3"), 
                             c("Oct 1 - Dec 31" = "4"))) %>% 
  group_by(year, season) %>% 
  summarise(sst = mean(sst, na.rm = T),
            sst_anom = mean(sst_anom, na.rm = T), 
            .groups = "keep") 

# Plot
quarter_summary %>% 
  ggplot(aes(year, sst_anom)) +
  geom_line(group = 1) +
  geom_point() +
  geom_smooth(method = "lm", 
              aes(color = "Regional Trend"),
              formula = y ~ x, se = F, linetype = 1) +
  stat_poly_eq(formula = y ~ x,
               color = gmri_cols("gmri blue"),
               aes(label = paste(..eq.label.., ..rr.label.., sep = "~~~")),
               parse = T) +
  scale_color_manual(values = c("Regional Trend" = as.character(gmri_cols("orange")))) +
  labs(x = "", 
       y = expression("Sea Surface Temperature Anomaly"~degree~C),
       caption = "Regression coefficients reflect annual change in sea surface temperature.") +
  theme(legend.title = element_blank(),
        legend.position = c(0.875, 0.05),
        legend.background = element_rect(fill = "transparent"),
        legend.key = element_rect(fill = "transparent", color = "transparent")) +
  facet_wrap(~season)

Marine Heatwaves

For the figures below heatwave events were determined using the methods of Hobday et al. 2016 and implemented using the R package {heatwaveR}.

A marine heatwave is defined as a situation when seawater temperatures exceeds a seasonally-varying threshold (usually the 90th percentile) for at least 5 consecutive days. Successive heatwaves with gaps of 2 days or less are considered part of the same event. The heatwave threshold used below was 90%. The heatwave history for Gulf Of Maine is displayed below:

# Use function to process heatwave data for plotting
region_heatwaves <- pull_heatwave_events(region_timeseries, threshold = 90)

Interactive SST Timeline

For anything we wish to host on the web there is an option to display tables and graphs that are interactive. Interactivity allows users to pan, zoom, and highlight discrete observations.

# Grab data from the most recent year through present day to plot
last_year <- Sys.Date() - 365
last_year <- last_year - yday(last_year) + 1
last_yr_heatwaves <- region_heatwaves %>% 
  filter(time >= last_year)


# Get number of heatwave events and total heatwave days for last year
# How many heatwave events:
last_full_yr <- last_yr_heatwaves %>% filter(year(time) == year(last_year))
num_events   <- max(last_full_yr$mhw_event_no, na.rm = T) - min(last_full_yr$mhw_event_no, na.rm = T)

# How many heatwave days
num_hw_days <- sum(last_full_yr$mhw_event, na.rm = T)
# Plot the interactive timeseries
last_yr_heatwaves %>% 
  filter(time >= last_year) %>% 
  plotly_mhw_plots()

Static SST Timeline

# Set colors by name
color_vals <- c(
  "Sea Surface Temperature" = "royalblue",
  "Heatwave Event"          = "darkred",
  "Cold Spell Event"        = "lightblue",
  "MHW Threshold"           = "coral3",
  #"MHW Threshold"           = "gray30",
  "MCS Threshold"           = "skyblue",
  #"MCS Threshold"           = "gray30",
  "Daily Climatology"        = "gray30")


# Set the label with degree symbol
ylab <- expression("Sea Surface Temperature"~degree~C)



# Plot the last 365 days
ggplot(last_yr_heatwaves, aes(x = time)) +
    geom_segment(aes(x = time, xend = time, y = seas, yend = sst), 
                 color = "royalblue", alpha = 0.25) +
    geom_segment(aes(x = time, xend = time, y = mhw_thresh, yend = hwe), 
                 color = "darkred", alpha = 0.25) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    scale_color_manual(values = color_vals) +
    scale_x_date(date_labels = "%b", date_breaks = "1 month") +
    theme(legend.title = element_blank(),
          legend.position = "top") +
    labs(x = "", 
         y = ylab, 
         caption = paste0("Climate reference period :  1982-2011"))

Heatwave Events

# Prep the legend title
guide_lab <- expression("Sea Surface Temperature Anomaly"~degree~C)

# Set new axis dimensions, y = year, x = day within year
# use a flate_date so that they don't stair step
base_date <- as.Date("2000-01-01")
grid_data <- region_heatwaves %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date))


# Set palette limits to center it on 0 with scale_fill_distiller
limit <- max(abs(grid_data$sst_anom)) *c(-1,1)


# Assemble heatmap plot
heatwave_heatmap <- ggplot(grid_data, aes(x = flat_date, y = year)) +
  
  # background box fill for missing dates
  geom_rect(xmin = base_date, xmax = base_date + 365, 
            ymin = min(grid_data$year) - .5, ymax = max(grid_data$year) + .5, 
            fill = "gray75", color = "transparent") +
  
  # tile for sst colors
  geom_tile(aes(fill = sst_anom)) +
  
  # points for heatwave events
  geom_point(data = filter(grid_data, mhw_event == TRUE),
             aes(x = flat_date, y = year), size = .25)  +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0)) +
  scale_y_continuous(limits = c(1980.5, 2021.5), expand = c(0,0)) +
  labs(x = "", 
       y = "",
       "\nClimate reference period : 1982-2011") +
  
  #scale_fill_gradient2(low = "blue", mid = "white", high = "red") +
  scale_fill_distiller(palette = "RdYlBu", na.value = "transparent", 
                       limit = limit) +
  
  #5 inches is default rmarkdown height for barheight
  guides("fill" = guide_colorbar(title = guide_lab, 
                                 title.position = "right", 
                                 title.hjust = 0.5,
                                 barheight = unit(4.8, "inches"), 
                                 frame.colour = "black", 
                                 ticks.colour = "black")) +  
  theme_classic() +
  theme(legend.title = element_text(angle = 90))





# Assemble pieces
heatwave_heatmap

SST Anomaly Maps

The 2021 global sea surface temperature anomalies have been loaded and displayed below to visualize how different areas of the ocean experience swings in temperature.

# Access information to netcdf on box
nc_year <- "2021"
anom_path <- str_c(oisst_path, "annual_anomalies/1982to2011_climatology/daily_anoms_", nc_year, ".nc")

# Load 2020 as stack
anoms_2020 <- stack(anom_path)

Heatwave Progression

Looking specifically at the last heatwave event, we can step through how the event progressed over time, and developing pockets or warmer/colder water masses.

# Identify the last heatwave event that happened
last_event <- max(region_heatwaves$mhw_event_no, na.rm = T)

# Pull the dates of the most recent heatwave
last_event_dates <- region_heatwaves %>% 
  filter(mhw_event_no == last_event) %>% 
  pull(time)


# Buffer the dates, start 7 days before
event_start <- (min(last_event_dates) - 7)
event_stop  <- max(last_event_dates)
date_seq <- seq.Date(from = event_start,
                     to   = event_stop,
                     by   = 1)


# Expand the area out to see the larger patterns
crop_x <- crop_x + c(-2.5, 2.5)
crop_y <- crop_y + c(-1.5, 1)

# Load the heatwave dates
data_window <- data.frame(
  time = c(min(date_seq) , max(date_seq) ),
  lon  = crop_x,
  lat  = crop_y)


# Pull data off box
hw_stack <- oisst_window_load(data_window = data_window, 
                              anomalies = T)


#drop any empty years that bug in
hw_stack <- hw_stack[map(hw_stack, class) != "character"]


##### Format the layers and loop through the maps  ####

# Grab only current year, format dates
this_yr   <- stack(hw_stack)
day_count <- length(names(this_yr))
day_labs  <- str_replace_all(names(this_yr), "[.]","-")
day_labs  <- str_replace_all(day_labs, "X", "")
day_count <- c(1:day_count) %>% setNames(day_labs)

# Progress through daily timeline to indicate heatwave status and severity
hw_timeline <- region_heatwaves %>% 
  filter(time %in% as.Date(day_labs))
####  Plot Settings:

# Set palette limits to center it on 0 with scale_fill_distiller
limit <- c(max(values(this_yr), na.rm = T) * -1, 
           max(values(this_yr), na.rm = T) )


# Plot Heatwave 1 day at a time as a GIF
day_plots <- imap(day_count, function(date_index, date_label) {
  
  # grab dates
  heatwaves_st  <- st_as_stars(this_yr[[date_index]])
  
  #### 1. Map the Anomalies in Space
  day_plot <- ggplot() +
    geom_stars(data = heatwaves_st) +
    geom_sf(data = new_england, fill = "gray90", size = .25) +
    geom_sf(data = canada, fill = "gray90", size = .25) +
    geom_sf(data = greenland, fill = "gray90", size = .25) +
    geom_sf(data = region_extent, 
            color = gmri_cols("gmri blue"), 
            linetype = 2, size = 1,
            fill = "transparent") +
    scale_fill_distiller(palette = "RdYlBu", 
                         na.value = "transparent", 
                         limit = limit) +
    map_theme +
    coord_sf(xlim = crop_x, 
             ylim = crop_y, 
             expand = T) +
    guides("fill" = guide_colorbar(
      title = expression("Sea Surface Temperature Anomaly"~degree~C),
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3, "in"),
      frame.colour = "black",
      ticks.colour = "black")) 
  
  
  
  
  #### 2.  Plot the day and the overall anomaly to track dates
  date_timeline <- ggplot(data = hw_timeline, aes(x = time)) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    scale_color_manual(values = color_vals) +
    
    # Animated Point /  line
    geom_point(
      data = filter(hw_timeline, time == as.Date(date_label)),
      aes(time, sst, shape = factor(mhw_event)), 
      color = gmri_cols("gmri blue"), 
      size = 3, show.legend = FALSE) + 
    geom_vline(data = filter(hw_timeline, time == as.Date(date_label)),
               aes(xintercept = time), 
               color = "gray50", 
               size = 0.5,
               linetype = 3,
               alpha = 0.8) + 
    labs(x = "", 
         y = "",
         color = "",
         subtitle = expression("Regional Temperature"~degree~"C"),
         shape = "Heatwave Event") +
    theme(legend.position = "bottom")
  
  
  ####  3. Assemble plot(s)
  p_layout <- c(
    area(t = 1, l = 1, b = 2, r = 8),
    area(t = 3, l = 1, b = 8, r = 8))
  
  # plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3))
  plot_agg <- date_timeline + day_plot + plot_layout(design = p_layout)
  
  
  return(plot_agg )
  
  
})


walk(day_plots, print)

Front Progression

Same idea as above but looking at the Belkin O’Reilly fronts rather than absolute values.

# Use belkin fronts function to get the sst fronts
this_yr_fronts <- map(unstack(this_yr), get_belkin_fronts) %>% 
  stack() %>% 
  setNames(names(this_yr))

# Set palette limits to center it on 0 with scale_fill_distiller
limit <- c(max(values(this_yr_fronts), na.rm = T) * -1, 
           max(values(this_yr_fronts), na.rm = T) )


# Build Plots for Animation

# Plot Heatwave 1 day at a time as a GIF
front_plots <- imap(day_count, function(date_index, date_label) {
  
  # grab dates
  sst_fronts_st  <- st_as_stars(this_yr_fronts[[date_index]])
  
  #### 1. Map the Anomalies in Space
  day_plot <- ggplot() +
    geom_stars(data = sst_fronts_st) +
    geom_sf(data = new_england, fill = "gray90", size = .25) +
    geom_sf(data = canada, fill = "gray90", size = .25) +
    geom_sf(data = greenland, fill = "gray90", size = .25) +
    geom_sf(data = region_extent, 
            color = gmri_cols("gmri blue"), 
            linetype = 2, size = 1,
            fill = "transparent") +
    scale_fill_distiller(palette = "RdYlBu", 
                         na.value = "transparent",
                         limit = limit) +
    map_theme +
    coord_sf(xlim = crop_x, 
             ylim = crop_y, 
             expand = T) +
    guides("fill" = guide_colorbar(
      title = "Front Strength?",
      title.position = "top",
      title.hjust = 0.5,
      barwidth = unit(3, "in"),
      frame.colour = "black",
      ticks.colour = "black")) 
  
  
  
  
  #### 2.  Plot the day and the overall anomaly to track dates
  date_timeline <- ggplot(data = hw_timeline, aes(x = time)) +
    geom_line(aes(y = sst, color = "Sea Surface Temperature")) +
    geom_line(aes(y = hwe, color = "Heatwave Event")) +
    geom_line(aes(y = cse, color = "Cold Spell Event")) +
    geom_line(aes(y = mhw_thresh, color = "MHW Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = mcs_thresh, color = "MCS Threshold"), lty = 3, size = .5) +
    geom_line(aes(y = seas, color = "Daily Climatology"), lty = 2, size = .5) +
    scale_color_manual(values = color_vals) +
    
    # Animated Point /  line
    geom_point(
      data = filter(hw_timeline, time == as.Date(date_label)),
      aes(time, sst, shape = factor(mhw_event)), 
      color = gmri_cols("gmri blue"), 
      size = 3, show.legend = FALSE) + 
    geom_vline(data = filter(hw_timeline, time == as.Date(date_label)),
               aes(xintercept = time), 
               color = "gray50", 
               size = 0.5,
               linetype = 3,
               alpha = 0.8) + 
    labs(x = "", 
         y = "",
         color = "",
         subtitle = expression("Regional Temperature"~degree~"C"),
         shape = "Heatwave Event") +
    theme(legend.position = "bottom")
  
  
  ####  3. Assemble plot(s)
  p_layout <- c(
    area(t = 1, l = 1, b = 2, r = 8),
    area(t = 3, l = 1, b = 8, r = 8))
  
  # plot_agg <- (date_timeline / day_plot) + plot_layout(heights = c(1, 3))
  plot_agg <- date_timeline + day_plot + plot_layout(design = p_layout)
  
  
  return(plot_agg)
  
  
})


walk(front_plots, print)

Warming Rate Ranking

If we look at the rates of change for each grid cell, rather than the observed temperature, its possible to rank how hot each location on earth is warming relative to others.

If we take the average ranking within the Gulf of Maine we can obtain the average warming rank for the area compared to the rest of the globe.

# 1. Warming Rates and Rankings
rates_path <- paste0(oisst_path, "warming_rates/annual_warming_rates")
rates_stack_all <- stack(str_c(rates_path, "1982to2020.nc"), 
                         varname = "annual_warming_rate")
ranks_stack_all <- stack(str_c(rates_path, "1982to2020.nc"), 
                         varname = "rate_percentile")

# Clip and average the warming rates for each area
mask_ranks <- function(mask_shape){
  
  # Get stats from ranks
  m1 <- mask(rotate(ranks_stack_all), mask_shape)
  m1 <- crop(m1, mask_shape)
  rank_mean <- cellStats(m1, mean)
  rank_min <- cellStats(m1, min)
  rank_max <- cellStats(m1, max)
  
  # Get stats from rates
  m2 <- mask(rotate(rates_stack_all), mask_shape)
  m2 <- crop(m2, mask_shape)
  rate_mean <- cellStats(m2, mean)
  rate_min <- cellStats(m2, min)
  rate_max <- cellStats(m2, max)
  
  # put in table
  table_out <- tibble("Mean Rank" = rank_mean,
                "Min Rank"  = rank_min,
                "Max Rank"  = rank_max,
                "Mean Rate" = rate_mean,
                "Min Rate"  = rate_min,
                "Max Rate"  = rate_max) %>% 
    mutate_all(round, 3)
  
  # spit them out
  return(table_out)
}


# Get the rank information that go with each area
region_ranks <- mask_ranks(region_extent)

# Prep it for text input.
avg_rank <- region_ranks$`Mean Rank` *100
avg_rate <- region_ranks$`Mean Rate`
low_rank <- region_ranks$`Min Rank` *100
low_rate <- region_ranks$`Min Rate`
top_rank <- region_ranks$`Max Rank` *100
top_rank <- ifelse(top_rank == 100, "greater than or equal to 99.5", top_rank)
top_rate <- region_ranks$`Max Rate`

Based on data from 1982-2020, the warming rates of Gulf Of Maine have been some of the highest in the world. The area as a whole has been increasing at a rate of 0.041\(^{\circ}C/year\) which is faster than 94.3% of the world’s oceans.

Over that same period locations within the Gulf Of Maine have been warming at rates as low as 0.014\(^{\circ}C/year\) and as rapidly as 0.088\(^{\circ}C/year\), corresponding to ranks as low as 50.8% and as high as greater than or equal to 99.5%.

Mapped below are the corresponding warming rates and their global rankings.

# Full plot
rates_map <- ggplot() +
  geom_stars(data = st_as_stars(rotate(rates_stack_all))) +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = gmri_cols("gmri blue"), 
          fill = "transparent", linetype = 2, size = 0.5) +
  scale_fill_viridis_c(option = "plasma") +
  map_theme +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = T) +
  guides("fill" = guide_colorbar(
    title = expression("Annual Temperature Change"~~degree~C~"/year"),
    title.position = "top", 
    title.hjust = 0.5,
    barwidth = unit(2.5, "in"),
    frame.colour = "black",
    ticks.colour = "black"))

# ranks map
ranks_map <- ggplot() +
  geom_stars(data = st_as_stars(rotate(ranks_stack_all))) +
  geom_sf(data = new_england, fill = "gray90", size = .25) +
  geom_sf(data = canada, fill = "gray90", size = .25) +
  geom_sf(data = greenland, fill = "gray90", size = .25) +
  geom_sf(data = region_extent, 
          color = gmri_cols("gmri blue"), 
          fill = "transparent", linetype = 2, size = 0.5) +
  scale_fill_viridis_c(option = "plasma",
                       limit = c(0.9, 1), oob = scales::oob_squish) +
  map_theme +
  coord_sf(xlim = crop_x, 
           ylim = crop_y, expand = T) +
  guides("fill" = guide_colorbar(
    title = expression("Global Percentile of Warming Rates"),
    title.position = "top", 
    title.hjust = 0.5,
    barwidth = unit(2.5, "in"),
    frame.colour = "black",
    ticks.colour = "black")) +
  labs(caption = "Ranking color scale truncated to display ranges of 0.9-1. Lower values will display as 0.9 or 90%")

# plot both
rates_map | ranks_map

Shifting Baselines

In 2021 NOAA is transitioning standard climatologies from the 30-year period of 1982-2011 to a new period spanning 1992-2020. Changes in climate regimes often does not reult in a uniform upward or downward change that is consistent throughout the year. The plot below shows just how both the average temperature, as well as the annual highs and lows have shifted.

When looking specifically at Gulf Of Maine here is how the expected temperature for each day of the year has shifted:

# Run heatwave detection using new climate period
heatwaves_91 <- pull_heatwave_events(region_timeseries, 
                                     threshold = 90, 
                                     clim_ref_period = c("1991-01-01", "2020-12-31")) 



# Subtract old from the new
heatwaves_91 <- heatwaves_91 %>% 
  mutate(clim_shift = seas - region_heatwaves$seas,
         upper_shift = mhw_thresh - region_heatwaves$mhw_thresh,
         lower_shift = mcs_thresh - region_heatwaves$mcs_thresh)


# Plot the differences
heatwaves_91 %>% 
  filter(time >= last_year) %>% 
  mutate(year = year(time),
         yday = yday(time),
         flat_date = as.Date(yday-1, origin = base_date)) %>% 
  distinct(flat_date, .keep_all = T) %>%        
  ggplot(aes(x = flat_date)) +
  geom_line(aes(y = clim_shift, color = "Mean Temperature Shift")) + 
  geom_line(aes(y = upper_shift, color = "MHW Threshold Change")) + 
  geom_line(aes(y = lower_shift, color = "MCS Threshold Change")) + 
  labs(x = "", 
       y = expression("Shift in Expected Temperature"~degree~C),
       color = "") + 
  theme(legend.position = "bottom") +
  scale_color_gmri() +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", expand = c(0,0))

 

A work by Adam A. Kemberling

Akemberling@gmri.org